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Abstract: Starting from the radiative transport equation we derive the 
scaling relationships that enable a single Monte Carlo (MC) simulation 
to predict the spatially- and temporally-resolved reflectance from ho- 
mogeneous semi-infinite media with arbitrary scattering and absorption 
coefficients. This derivation shows that a rigorous application of this single 
Monte Carlo (sMC) approach requires the rescaling to be done individually 
for each photon biography. We examine the accuracy of the sMC method 
when processing simulations on an individual photon basis and also demon- 
strate the use of adaptive binning and interpolation using non-uniform 
rational B-splines (NURBS) to achieve order of magnitude reductions in 
the relative error as compared to the use of uniform binning and linear 
interpolation. This improved implementation for sMC simulation serves as 
a fast and accurate solver to address both forward and inverse problems and 
is available for use at http : / /www . virtualphotonics . org/. 
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1. Introduction 

Real-time optical analysis of turbid media requires rapid radiative transport solvers to predict 
and analyze remitted light signals for the non-invasive determination of optical scattering and 
absorption characteristics. While analytical techniques are often preferred for their simplicity 
and speed, they are often based on the diffusion approximation to the radiative transport equa- 
tion (RTE). The limitations of the diffusion approximation are well known and result in accu- 
rate radiative transport estimates only when the reduced scattering coefficient, jl f s = jl s (l — g), 
is much larger than the absorption coefficient \i a and over spatial scales that are larger than the 
transport mean free path, /* = (ji a + [1,2]. 

These limitations have motivated extensive efforts to solve the RTE, using either stochastic 
solvers such as Monte Carlo [3] or deterministic solvers via analytical [4] or numerical [5] ap- 
proaches. RTE solutions provide accurate radiative transport predictions regardless of optical 
properties and on spatial scales approaching a single scattering length; t s = \j\i s ~ 10— 50/im 
for biological tissues in the red and near-infrared spectral regions. Unfortunately, these RTE 
solution approaches are computationally intensive and currently impractical for real-time ap- 
plication. For the case of homogeneous turbid media, several groups have proposed methods 
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that utilize the results of a single MC simulation to provide rapid transport-rigorous RTE solu- 
tions for a class of related radiative transport problems. 

Graaff and co-workers first investigated the use of a single MC (sMC) simulation to obtain 
results for turbid slabs with different sets of optical properties [10]. They demonstrated that the 
trajectories of unabsorbed photons in a single Monte Carlo simulation can be used to estimate 
the total reflectance and transmittance of turbid slabs for various single-scattering albedos, a = 
Us I \H>a + Ms) by recording the number of photon interactions within the simulated medium. 
Using this approach they predicted the spatially-resolved reflectance for different albedo values 
in cases for which the interaction coefficient, jl t = fl a + jl s , was held constant. 

Kienle and Patterson (KP) extended the capabilities of sMC simulations to provide the 
spatially- and temporally-resolved reflectance R(r,t) for a homogeneous semi-infinite medium 
with arbitrary optical properties [11]. This sMC method is based on scaling and re- weighting 
a discrete representation of the time-resolved reflectance provided by a reference MC simula- 
tion in a non-absorbing semi-infinite medium. In this approach the evaluation of R(r,t) with 
arbitrary optical properties, requires the interpolation of the discrete representation of R(r,t) in 
both the spatial and temporal dimensions. While powerful, the discrete representation of the 
reference R(r,t) and interpolation both introduce errors that are often amplified when the scal- 
ing is applied. Most recently, Alerstam and co-workers improved upon the KP sMC method 
by applying the scaling relationships on a photon-by-photon basis prior to spatial and temporal 
binning [13]. This approach requires storage of the radial exit position and total path length of 
each detected photon, which are individually reprocessed to evaluate the photon weights for 
media with different optical properties. 

The motivation and objectives of this manuscript are two-fold. While sMC is widely used 
by the biomedical community [12, 14-17], a theoretical analysis that provides a clear line of 
derivation from the RTE to the scaling relations that form the basis of sMC has not been pre- 
sented. Thus it is unclear whether differences between results obtained from conventional MC 
vs. sMC methods arise from stochastic variations, shortcomings of the scaling relationships 
and/or inaccuracies associated with the binning and interpolation used in the sMC implemen- 
tation. To address these issues, our first goal is to determine to what extent the sMC scaling 
relations are transport-rigorous and to provide any ancillary conditions necessary for their va- 
lidity. Second, in order to establish 'best practices', we aim to examine the impact of binning 
and interpolation strategies on the relative error/variance of the sMC predictions. 

2. Theory 

The temporally- and spatially -resolved reflectance R(r,t) is sufficient to derive all reflectance 
metrics typically used to characterize a turbid system, i.e., steady-state, temporal frequency- 
domain, and/or spatial frequency-domain reflectance. While conventional Monte Carlo simu- 
lations can be performed without explicit discretization of the phase space, their interpretation 
usually requires a mesh to be imposed on the phase space. Thus, the detected photons are placed 
into finite bins to visualize the simulation output. For a homogeneous half space occupying 
z > 0 and possessing cylindrical symmetry, the spatially- and temporally-resolved reflectance, 
R(r,t) \ z= q- at a given location (ro,*o) on the tissue surface can be estimated by 

* (ro '' o) ~ AAAt i 0 J ro y § 2- L(r ' Q '^L=o-^' Q ) dQ2 ^ drd ^ (!) 

where L(r, £2, t) is the radiance, Q the unit direction vector, h the outward pointing unit normal 
vector, AA = n (Yq + Ar) — nr\ the surface area from which the photons are collected, and S 2_ 
the solid angle interval corresponding to the upward-pointing hemisphere. Many photons will, 
in general, exit a given cylindrical ring within the same time interval. Thus, when these photons 
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are gathered to obtain the value of R(ro, to), a discretization error occurs because the quadrature 
interval has been represented by the nominal values (ro,*o)- This error vanishes as the number 
of simulated photons increases without bound and as the integration intervals — >> 0. 

Two scaling relationships are used in the sMC method to obtain spatially- and temporally- 
resolved reflectance estimates for a homogeneous half space. The first uses the reflectance of a 
'reference' medium R r (r,t) with ji a = 0 and \l s = /l^ r to determine the reflectance R(r,t) from 
a second non-absorbing medium with arbitrary scattering coefficient jl s [11]: 



R(r,t) 



^Rr(^r^t). (2) 



The second scaling relationship relates the reflectance from a non-absorbing 'reference' 
medium (/i^ r = 0) to a second medium with identical scattering but non-zero absorption ji a : 

R(r,t) =R r (r,t)exp(-n a ct) (3) 

where c is the velocity of light propagation in the medium. The analysis that follows will estab- 
lish a rigorous foundation for these two scaling relationships. 

We begin our analysis by considering the source-free RTE for a non-absorbing medium: 



1 dL (r,Q,f) 
c dt 



+ VL(r,Q,f) Q = -n s L(r,Cl,t) + ^ / /?(^' Q)L(r,Q',f) dft' 



(4) 



where L (r, Q, f) is the radiance at a given location r at at given time nn a given direction Q. c is 
the velocity of light propagation in the medium and is equivalent to the ratio between the speed 
of light in vacuum and the medium's refractive index (= co/n). jl s and \l a are the scattering and 
absorption coefficients, respectively, and p(Cl —> Q) is the single- scattering phase function, a 
probability density function for scattering from Q to Q. For boundary conditions, we assume 
a vanishing radiance at remote locations in the medium 

L(r->oo ? Q ? f) ->0, (5) 

and allow a possible refractive-index mismatch at the surface: 

L(r, & r ,t)=& F (h- Q/)L(r, Q /5 1) , r e 5, h -Cli> 0, (6) 

where ^ is the angular distribution of the Fresnel reflectance for unpolarized light and Q r is 
the unit vector in the direction of the reflected ray for an incident direction Q/. Thus &p(n • Q/) 
represents the probability for internal reflection for a photon incident on an interface with a 
refractive index discontinuity. 

We first express the RTE and boundary conditions in dimensionless form using the following 
variables: 

/ T = fl s ct 

p = Us* (7) 
[A(pAr)=L(rAf) /L 0 , 

where T is a dimensionless time, p a dimensionless position and A (p,Q,r) a dimensionless 
radiance defined as the ratio between the radiance, L (r, Q,f) , and a characteristic radiance Lq. 
Substitution of these variables into the RTE yields: 



— ^ 7 Q? T ^ + VA (p,Q,r) • Q = -A (p , Q, t) + / p(o! ^Cl)A(p,Cl\T)dn'. 

uX «/S 2 



(8) 
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The two boundary conditions, Eqs. (5) and (6) become 

A(p ^oo,Q, T ) ^0, and (9) 

A(p A, t) = £% F [n • Q/)A(p , Q/, t) , peS, h Cli> 0. (10) 

Note that the non-dimensionalization does not alter the form of the boundary conditions. Hav- 
ing defined the boundary conditions, the diffuse reflectance is given by 

R(r,t) = j^_L(rAt)\ z=0 -{fi-Cl) dCl. (11) 

By analogy, a dimensionless reflectance is defined as 

R(p,T)= [ A(p,Q,t) (n-Q) dCl. (12) 

Js 2 ~ z=0 

To derive the first sMC scaling relation, Eq. (2), we consider the dimensionless analog of 
Eq.(l) 

R(po,T 0 ) = - J — / / A(p,Q,r) (n-O) dQ 2^p dp dr (13) 

where 

AA = 7r(po + Ap) 2 - Trpo 2 . (14) 

Since the ratio between the radiance L(r,Q,^) and the dimensionless radiance A(p,Q,r) is 
defined as Lo, one can easily verify that the ratio between R(r,t) and 7?(p,r) is also equal to 
Lo. The triple integral in Eq. (1) represents a characteristic energy Eo that is emitted in the 
space-time bin considered. This characteristic energy Eo is equal to the product of the number 
photons N emitted in the bin multiplied by the energy of each photon, hv. Because we integrate 
L(r, Q, t) and A(p, Q, t) over equivalent ranges of solid angle, space, and time, we can write 

rt 0 +At /T 0 +Ar r ^ />T 0 +At rp 0 +Ap r _ 

/ / / L(r,£l,t)d£l2nrdrdt = E 0 / / A(p,£2,r)dQ27rpdpdT. 

Jt 0 Jr 0 J$ 2 ~ Jt 0 Jpo JS 2 ' 

(15) 

Using this equivalence, as well as the definitions for R(ro,to) and 7?(po,To) in Eq. (1) and 
Eq. (13), we find that Lo is the ratio between the reflectance and the dimensionless reflectance 

tf(po,T 0 ) AAA? 

With this definition for Lo, the reflectance can be written in terms of the dimensionless re- 
flectance as 

R(r,t) = cn*EoR(p,i;) = cn*Eo [ A(p,Q,r) (h-Cl) dCl. (17) 

This equation shows that the diffuse reflectance for a specific value of the scattering coefficient 
is uniquely related to the dimensionless reflectance multiplied by a reference radiance Lo. 

Thus Eq. (2) can be formulated by considering a reference reflectance R r {r,t) for a non- 
absorbing medium with scattering coefficient jl s = /l^ r . Using Eq. (17), to express the reference 
reflectance within the KP sMC method in terms of the dimensionless reflectance, we get: 

R r (r,t) =cnl r E 0 R r (p,T) =cnl r EoR r (n s , r r,cn Sir t), (18) 
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while the diffuse reflectance for an arbitrary value of n s becomes 

R(r,t) = ctfW(p , T) = c^E 0 R I ^p r , T r ) ( 1 9) 

V l^s, r l^s, r J 



where p r = /l^ r r and % r = \i Sl rCt. By rewriting the dimensionless reflectance in the last equation 

as 

~( ll s ll s \ (ll s \ 1 n ~ 

R\ p n T r = R r \ r, 1 — v — (20) 

and substituting back into Eq. (19), we obtain Eq. (2) 




R(r,t) = -g-R r [-?-r,-?-t\, (21) 



which demonstrates the consistency of the first sMC scaling relationship with the RTE. 

The second scaling relation used within the sMC method, Eq. (3), relates the reflectance from 
a non-absorbing reference medium with the reflectance from a second medium possessing iden- 
tical scattering properties but finite absorption. To evaluate the impact of absorption, consider 
that every photon emitted from the target surface is characterized by a photon exit position r/, 
time of flight t[, and weight w/. While photon absorption is an analog phenomena, an effec- 
tive technique to reduce the variance of the resulting photon weights is to diminish the photon 
weight continuously along its trajectory according to Beer's law i.e., exp(— fl a cti) [21]. This 
approach, known as continuous absorption weighting, provides a rigorous unbiased stochastic 
estimator of the radiative transport process [3]. Consider a reference simulation with zero ab- 
sorption and calculate N r (ro,to) as the sum of the weights of the photons emitted from the top 
surface within a specific radial and temporal bin (ro,*o) 

Afr(r 0 ^o)=X w '' = X w o (22) 

where wo = 1 — R sp is the initial weight of the photon reduced by the specular reflection given 
by the refractive index mismatch between air and tissue. Thus, following Eq. (1), we can express 
the reference reflectance for the bin represented by ro and to as 

1 /*o+Af rro+Ar r ^ 

*'( r o>'o) = 777- / / / hvcJK r (r 0 ,to,n)\ z=0 - d&lnr&r&t, (23) 

where ^(ro,^?^) [ m_3 sr_1 ] represents the volumetric photon density per unit solid angle at 
location ro and time to in the reference MC simulation. 

For a MC simulation with \i a ^ 0 the weight of each photon must be decreased exponentially 
according to Beer's law, i.e., Af r (ro,*o) = = X^oexp (—f± a cto). Thus the reflectance is 

I rtQ+At rro+Ar r . . 

^fa'to) ^ 7777 / / /. hvc^K(ro,to,n)\ z=0 -Qxp(-ii a cto)(h-a) dQ2nrdrdt 
= R r (r 0j t 0 )exp(-ii a cto). 

(24) 

If we consider continuous spatial and temporal variables we obtain the expression 

R(r,t) =R r (r,t)txp(-ii a ct) (25) 

which recovers the second sMC scaling relationship Eq. (3). Equation (25) is rigorous only 
when applied individually to each emitted photon. That is, the value of t used must correspond 
to the photon's time-of-flight rather than to a nominal time of the temporal bin in which it is 
tallied. 
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3. Materials and Methods 



The flowcharts in Fig. 1 depict three approaches used here to obtain and compare MC estimates 
for spatially- and temporally-resolved reflectance from a homogeneous semi-infinite medium. 
In the first approach, a conventional Monte Carlo (cMC) simulation is used to provide "gold 
standard" reflectance estimates. The second method is based on the single Monte Carlo ap- 
proach applied on a photon-by-photon basis (sMC p ). The third approach is similar to the KP 
method where the reflectance for any value of the optical scattering and absorption, R(r,t), is 
obtained by scaling a reference signal, R r (r,t), obtained from a MC simulation with \i s = /i^ r 
and \± a = 0. As this latter method employs interpolation of the reference reflectance prior to 
application of the scaling equations, we denote it as sMQ. Below, we describe the detailed 
procedure used in each of these methods and also describe alternate binning and interpolation 
schemes that we have examined to potentially improve upon the KP sMQ method. 

3.1. Monte Carlo Simulations 

We modified the well-known MCML MC simulation package developed by Jacques and 
Wang [21] to provide time-resolved output. All simulations run for this study inject photons in 
a semi-infinite homogeneous medium with n = 1.4 using a collimated 'pencil' beam oriented 
perpendicular to the top surface. The Henyey-Greenstein phase function with an anisotropy co- 
efficient of g = 0.8 is used to model the angular distribution of single- scattering events. When 
a photon encounters the z = 0 surface, we select a pseudorandom number uniformly distributed 
in [0,1) and compare it to the Fresnel reflectance to determine whether the photon should exit 
or be internally reflected. To control the computation time, simulated photons are terminated 
before exiting the medium if they experience 3 x 10 5 collisions. To achieve good statistics, the 
propagation of 10 8 photons was simulated for each set of optical properties. 

3.1.1. Conventional Monte Carlo (cMC) 

In cMC simulations, the step size for photon propagation is sampled using the probability den- 
sity function s = — ln(^)/ji s , where jl s is the scattering coefficient and § is a pseudorandom 
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Fig. 1. Flow chart representing the processes used to obtain the spatially and temporally 
resolved reflectance from the different Monte Carlo approaches examined in this study. 
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number uniformly distributed in the interval [0,1). For each photon simulated, the exit posi- 
tion Y[ and total time of flight U are recorded. The weight of each photon w; is then calculated 
using continuous absorption weighting according to Beer's law for a specific value of jl a . To 
reduce the computational burden we terminated all photons that traveled to a depth exceeding 
1 m. After storing exit position, time, and weight (r/,f/, w/) of each of the reflected photons in a 
database, the recorded photon weights were averaged over a set of spatial and temporal bins to 
obtain the reflectance signal [13]. 

3.1.2. Single Monte Carlo on a photon-by-photon basis (sMC p ) 

The application of sMC on a photon-by-photon basis requires simulation of the photon trajecto- 
ries and storage of the photon exit locations and times in a database. We perform a dimension- 
less MC simulation with zero absorption, set jl s equal to a unitless value of one and sample the 
photon step size from s = — ln(<^ ) . The maximum depth for photon propagation was adjusted to 
match the termination condition of the cMC simulations. The dimensionless position and time 
(p/,T/) of each exiting photon was recorded. Since \l a = 0 in this simulation, all the emitted 
photons had the same weight w ; = 1 — R sp - To obtain the temporally- and spatially -resolved 
reflectance for any value of \l a and \l s we first apply Eq. (7) to determine the exit position and 
time, Vi and t{, and apply the scaling relationships, Eqs. (2) and (3) to each exiting photon to 
determine the proper photon weight. Once the exiting position, time, and weight of each photon 
is obtained, the individual photon weights are averaged over a set of spatial and temporal bins 
to obtain the reflectance signal as in the cMC simulations. 

3.1.3. Interpolated single Monte Carlo (sMQ) 

sMQ refers to the method established by Kienle and Patterson which uses a cMC simulation 
to form the reference reflectance R r {r,i) using reference optical properties j±' s r = 1 mm -1 and 
jl a = 0. The resulting photon database is then processed to obtain the reference reflectance for 
a set of nominal positions, R r {r^ti). The reference values sampled through the MC simulation 
are smoothed by fitting to an analytical function and then interpolated. The reflectance for any 
combination of optical properties is obtained by applying the scaling relations Eqs. (2) and (3) 
to the interpolated representation of the reflectance. 

3.2. Error Analysis 

Our objective is to compare both sMC p and sMQ relative to the conventional MC simulations. 
The relative error for the spatially- and temporally- resolved reflectance is given by: 

£ ^ {rk > tl) = ' (26) 

The cMC values chosen as gold standard are stochastic estimates of the RTE solution. The 
relative error estimates are meaningful only if the relative uncertainty associated with the 
R C Mc( r k,ti) estimates is at least an order of magnitude smaller. This relative uncertainty is 
calculated as 

_ / ,x GcMc(rk,ti) 

are ' (r ^ /) = w^y- (27) 

3.3. Binning 

To represent the MC data, the detected photon weights must be tallied in a set of radial and 
temporal bins of finite width. For each reflected photon emitted at radial position r ; and time 
ti, a tally of w; is added to the reflectance R{r^ti) where k is the index associated with the 
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radial ring that contains r; and / is the index for the time interval that includes t[. This process 
is necessary to obtain R(r^t\) using either the cMC or sMC p approaches. For the KP sMQ 
approach, the reference MC simulation is first binned to obtain the nominal values R{r^t\). 
The binned reference MC simulation is then processed to obtain R r (r,t) as illustrated in Fig. 1. 

Three different discretization techniques have been examined in this work. The first is com- 
monly known as equal width discretization (EWD) [25] and is based on the subdivision of the 
physical range into m r equispaced radial bins and m t equispaced time bins resulting in a total 
of m r x m t bins. This binning strategy is most commonly adopted for the evaluation of the re- 
flectance signal from the output of a cMC simulation [17]. A slight variation presented in the 
literature utilizes multiple sets of uniformly-sized bins [11]. 

The second approach, known as equal frequency discretization (EFD) [26], divides the sorted 
values into m r x m t intervals of irregular width so that each interval contains approximately 
the same number of samples. Instead of fixing the boundary of the bins a priori, the bins are 
built dynamically based on the radial and temporal cumulative distribution functions of the 
collected photons. Thus each interval contains N out / (m r x m t ) photons where N out is the number 
of photons that exited the phase space and were not terminated during the simulation. This 
binning process requires sorting of the photon database according to the radial position. Once 
the data set is split into m r bins in space with the same number of photons in each bin, the 
photons in each spatial bin are sorted according to the photon time of flight. The reflectance 
values are then computed by using EFD along the time dimension. Both EFD and EWD were 
used to provide a discrete representation of the reflectance R(rk,ti) from the cMC simulation 
output. For all the combinations of optical properties under investigation we save the bin limits 
obtained using EFD and use them to obtain Rfa^ti) from the photon database generated with 
sMC p . 

We developed a third binning strategy based on an adaptive width discretization (AWD) 
algorithm to construct the reference reflectance used for the sMQ. In AWD the number of 
bins is not imposed a priori and the bin limits are dynamically built to ensure that each bin 
has a sufficient sample size of photons to reduce variance while also spacing the bins so that 
they accurately capture any rapid spatial and/or temporal dynamics of the reflectance signal. To 
accomplish this, we established the following binning criteria: (/) a minimum ofrc r = 5x 10 5 
photons must occupy each radial bin; (ii) the number of photons occupying each time bin, 
n t (r), increases from 2500 to 5000 over the radial distance range; (Hi) the interval between two 
adjacent radial bins cannot exceed 0.2 mm i.e., (r^i — r^) < Ar max = 0.2 mm; (iv) the interval 
between two adjacent time bins cannot exceed 20 ps, i.e., (ti + \ —r{)< At max = 20 ps; and (v) 
the relative change in the reflectance value between two adjacent time bins cannot exceed 10% 
i.e., (Rk+i —R^/Rjc < AR re i max 0.1. While other values for n r ,n t (r), Ar max ,At max and AR min 
were tested, the values above provided a smaller variance for the reflectance values R(rk,ti) and 
reduced noise in the data while accurately capturing sharp features of the reflectance signal. 

Independent from any specific binning algorithm, the representation of all the photons col- 
lected within a bin with nominal coordinates (r^ti) introduces a discretization error. Typically 
each bin is represented by its midpoint. To reduce this error we used the mean of the radial po- 
sition and time of flight of all the photons collected within each bin. While this error reduction 
strategy has not been examined carefully, we noted that when large bins are present in EFD or 
EWD, the use of the midpoint leads to artifacts in the measured reflectance. Regardless of the 
binning approach used, the value of the reflectance calculated over an annular bin of area AA^ 
and time span Ati is computed as [21] 

R ^=ikk X N% Wi > (28) 
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where N is the total number of injected photons. The weight w; is set to zero for all the pho- 
tons collected outside the bin under consideration. The standard deviation of the reflectance 
measured in the specific bin is calculated as [27]: 



1 



N 




1/2 



I 



-R 2 (r h t t ) 



(29) 



i=l 



3.4. Reference Reflectance 

In principle, the use of Eqs. (2) and (3) enables a continuous representation of R(r,t) for any 
absorption and scattering coefficient. However when utilizing the results of a cMC simulation 
with a finite number of photons, it is possible to measure the reflectance only for a discrete set 
of positions. Thus, it is necessary to interpolate between these points to provide a continuous 
representation of R(r,t) for an arbitrary set of optical properties. To represent the reflectance 
signal over the physical domain of interest, we establish a native set of reflectance values, 
Rr(ric,ti) over a time range from 0 to 20 ns and a spatial range from 0 to 100 mm. 

The adaptive binning technique described above is used to provide a representation of the 
reference simulation output which must subsequently be interpolated. Formation of an interpo- 
lating surface requires a support mesh of points. Since the number of the reflectance values and 
their temporal locations are different for every radial bin, we resample each reflectance curve 
R(rk = constant,^) over a uniform time vector T to obtain a regular grid that is more suitable 
for interpolation. These time points are chosen in a way that reflects the cumulative distribution 
function of all the photons, using a decreasing density of sampled points per unit time in the 
range [0-20] ns. 

3.5. Interpolation and Resampling 

The AWD approach had one limitation in that small numbers of photons were typically col- 
lected at late times. For this reason the last time coordinate value that satisfied condition (//), 
R{ r kitiast)i occurs at a time smaller than 20ns for radial distances r# < 15 mm. To obtain nom- 
inal values for the reflectance in this time range we defined a set of bins 100 ps in width from 
this valid time location to 20 ns. For each radial bin that exhibited this problem, the photons 
collected within its limits are tallied in these bins. The variance associated with the photon col- 
lected in these 100-ps- width bins is still too large to use them directly for the creation of the 
reference surface. Thus we fit these points to the following nonlinear function, f t (t), with three 
free parameters (a,b,c) 



This function has the form of a source-image Green's function without absorption decay within 
the standard diffusion approximation [2]. This functional form fits the data very well most 
likely because the photons collected at these remote distances and long times have undergone 
a very large number of collisions. Figure 2(a) shows the result of this approach to represent the 
reference reflectance for r& = 1 mm. 

To improve the accuracy in the representation of the initial rise of the reference reflectance 
signal we introduce additional time points at earlier times via extrapolation. For every radial bin 
we record the arrival time of the first photon to. The time span from to to the nominal position 
of the first bin, t\ , was split into ten uniform bins in which the photons were tallied to provide 
additional reflectance values. The logarithm of the reflectance values sampled with this method 
is fitted to a third degree polynomial. For each radial bin, this polynomial is used to extrapolate 
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Fig. 2. (a) Results of the fitting approach used to extrapolate points on the tail of the time- 
resolved reflectance R(t) for = 1mm. (b) Values of R r {r^t{) for ry- e [5 — 100] mm 
obtained after binning, extrapolating and resampling. 



values for the points of T between the arrival time of the first photon, to, and the time of the 
first reflectance value obtained with the binning process, t\ . 

To resample the reflectance curves at the temporal locations specified by T, we fit a 3rd 
degree Non Uniform Rational B-spline (NURBS) [28,29] curve to each 7?(r^ = const, f/) curve. 
After calculating the approximating NURBS curve we resample it to obtain R r (rk = const, t\ = 
T). The values of R r {r^t{) obtained after resampling is the final set used in the sMQ method 
and are shown in Fig. 2(b) for radial distances spanning 5 to 100 mm. The noise reduction 
obtained by resampling allowed the use of the reference points R r (rk,ti) as native points for the 
interpolation necessary to calculate R(r,t) for continuous values of r and t for an arbitrary set 
of optical properties using Eqs. (2) and (3). We use a 3rd degree NURBS interpolation in both 
position and time to represent R r {r^t{) as a 2D surface. 

4. Results 

4.1. cMC Simulations and Effects of Discretization 

We first consider the relative standard deviation of the reflectance provided by our "gold stan- 
dard" cMC simulations as these are the results against which the sMC approaches will be com- 
pared. We calculate the relative standard deviation using Eq. (27) after binning the reflected 
photons provided by the cMC simulation using either EFD or EWD. We used a set of 200x200 
bins to calculate the reflectance values over spatial and temporal ranges of r e [0, 50] mm and 
t G [0, 5] ns, respectively. We performed cMC simulations for 15 sets of optical properties where 
the absorption and reduced scattering coefficients were in the range \l a £ [10 _3 ,0.3] mm -1 and 
H f s £ [0.5, 2] mm" 1 

Figure 3(a) provides a colored contour plot of spatial and temporal dependence of the rela- 
tive standard deviation of the reflectance estimates o re i(r,t) for the case of \l' s — 1mm -1 and 
\± a = 0.1mm -1 . Figure 3 shows that the values of o re i are remarkably uniform over the en- 
tire domain. Across all 15 sets of optical properties examined, the mean of the relative standard 
deviation across all space and time using EFD was less than 3%. For larger values of the absorp- 
tion coefficient, a re \ increases for combinations of short distances and late times and for long 
distances and early times. This is because the reflectance signal is weaker in these regions and 
the span covered by the EFD bins is large at late times in order to collect the required number 
of photons. The relative standard deviation at any given position/time pair exceeds 5% only for 
jl a > 0.1 mm -1 . By contrast, in Fig. 3(b) we see that the relative standard deviation obtained 
with EWD is not uniform and increases strongly for longer times and larger distances with 
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Fig. 3. Relative standard deviation of the reflectance provided by discretizing the results of 
a cMC simulation for fi a = 0.1 mm -1 and ji f s = 1 mm -1 using (a) 200 x 200 EFD bins, (b) 
200 x 200 EWD bins, (c) 250 x 1000 EWD bins. 

a mean value larger than 5%. Moreover, for all the sets of optical properties examined, there 
were position/time pairs for which the relative standard deviation in the reflectance estimates 
exceeded 100%. 

A comparison of the uncertainty obtained using the same number of bins using the two 
techniques may lead one to underappreciate the improvement obtained with EFD. It is true that 
in certain regions, the uncertainty in the cMC estimates obtained when using EWD with a 200 
x 200 set of uniform bins is smaller than that obtained using EFD. This is because the uniform 
bins have radial and temporal width of Ar = 0.25 mm and At = 25 ps. Thus for combinations 
of small radial positions and times, the number of photons tallied within these limits is very 
large and the statistical noise is reduced. However bins that are 25 ps in width are too large 
to capture the initial peak of the reflectance signal. For a more fair comparison we use EWD 
with 250 radial intervals of width Ar = 0.2 mm and 1000 time bins of width At = 5 ps. These 
results are given in Fig. 3(c) and shows that the performance of the EWD at large times and 
radial positions is adversely effected when the number of bins is increased. In fact for this 
set of 250 x 1000 bins the mean a re \ across all the time and spatial bins exceeds 15% for all 
the sets of optical properties examined. Thus while EWD is simple to implement, the relative 
standard deviation will exhibit significant variability since the number of photons tallied in 
each of the bins varies widely. EFD, on the other hand, reduces the variation in the number 
of tallied photons across all bins by imposing finer sampling in the regions of phase space 
where number density of reflected photons is large while also imposing coarser discretization 
in regions where the reflected signal is sparse. This has the overall effect of not only reducing 
the relative standard deviation but also also minimizing its variation over the phase space. 

4.2. Comparison of single Monte Carlo Techniques 

Now with the relative standard deviation of the "gold- standard" cMC simulations characterized, 
we can determine the relative error provided by sMC p and sMQ approaches established in §3 
as compared to the reflectance provided by a cMC simulation R C Mc( r k^i) represented by a set 
of 200 x 200 radial/temporal bins formed using EFD. 
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For all the combinations of optical properties under investigation, we process the sMC p data 
set to evaluate the values (r ; ,^,w ; ) from (Pj,T;) using the scaling imposed by the value of fl s 
and modifying the photon weight in media with \i a ^ 0 by applying Beer's law using the values 
of the photon exit time t[. To calculate R S MC p ( r k^i) for the same nominal spatial and temporal 
locations, we use the bins defined by the cMC simulation for the specific combination of the 
optical properties to tally the photon weights and build up the reflectance values. We use the 
NURBS surface obtained by the procedure described in §3.4 to calculate R S MCi( r k^i)- First, 
the scaling given by the scattering coefficient is imposed using Eq. (2). Once established, the 
effects of non-zero absorption is calculated using Eq. (3). 

Figures 4, 5, and 6 show (a) the reflectance determined from the cMC simulations on a 
logarithmic scale along with the relative error in the reflectance predictions provided by the (b) 
sMCp and (c) sMQ approaches for different combinations of optical properties. The spatial and 
temporal ranges shown extend only to 20 mm and 1 ns to provide reasonable resolution over 
the region where the largest contribution to the reflectance signal resides. Figure 4(b) shows 
that the sMC p method has a relative error of < 3% over the entire spatial and temporal domain 
considered. As this is roughly equal to the maximum relative standard deviation of the 200 x 
200 EFD representation of the cMC, we can conclude that the reflectance estimates obtained 
using the cMC and the sMC^ methods are equivalent. In fact, for all sets of optical properties 
investigated, the mean value of the relative error £ r el ,sMC p ( r k,ti) over the entire physical range 
(r G [0,50] mm and t G [0,5] ns) never exceeds 3% which is smaller than the mean of the 
relative standard deviation associated with the cMC measurements. Thus there appears to be 
no statistically significant difference between the reflectance estimates given by cMC and sMC p 
simulations when using a set of 200 x 200 radial/temporal bins formed using EFD. 

The relative error of the reflectance calculated using the sMQ approach displays two char- 
acteristic features regardless of optical properties. Firstly, the relative error in the first temporal 
bin can be as large as 100%. It is important to note that this large error is measured only for the 
first time bin, a location where the reflectance value is much smaller then the peak value. Sec- 
ondly, the relative standard deviation associated with the first time bin at large source detector 
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Fig. 5. (a) Reflectance signal obtained with EFD from the cMC simulation, (b) relative 
error of the reflectance provided by the sMC^ method, (c) relative error of the reflectance 
provided by the sMQ method, for \i a = 0.01 mm -1 and ji f s = 1 mm -1 . 

separations is generally larger than the average error measured over the entire physical range. 

We also see specific features of the relative error depending on optical properties. For ex- 
ample, for ji a values smaller than 0.1 mm -1 across all jl' s values, the mean relative error of 
the sMQ approach tallied over the entire temporal and radial range excluding the first tempo- 
ral bin is smaller than that obtained with sMC p method. This reflects the advantage gained by 
interpolating the reference signal which reduces the impact of statistical noise. For \i a values 
equal to 0.1 mm -1 and greater, we see that while the mean relative error for the reflectance 
values provided by the sMQ approach is larger than that provided by the sMC^ approach, the 
error itself remains below 5%. From Figs. 4(c) and 5(c) one can see that the error introduced 
by using Eq. (30) to extrapolate the reflectance at long times is very small for small jl a values. 
The values obtained in the latest time bins for short distances are obtained by scaling the ref- 
erence reflectance in a range where the nominal reference values were obtained with the fitting 
approach, and the relative error remains very small. 

For higher \i a values, the error measured in the latest time bins at small radial locations 
increases as shown in Fig. 6(c). the error exceeds 50% for ji a > 0.1 mm -1 and approaches 
100% for \i a = 0.3 mm -1 over the physical range shown. This large error is due to the fact 
that Eq. (3) is approximate when applied to binned data where a median value of t is used to 
rescale the weights of all photons captured in a discrete time bin. It is important to note that the 
relative standard deviation of these bins is also large. While these relative error values are large, 
one must recognize that the reflectance values span a significantly larger dynamic range for \i a 
values of 0.1 mm -1 and larger. Specifically, within the radial and temporal range considered, 
the reflectance values span at least 8 orders of magnitude. When viewed in this context, we see 
that these large relative errors occur in regions where the reflectance is already very small. 

Since the reference signal represented by the NURBS surface extends to temporal and spa- 
tial values of 20 ns and 100 mm, respectively, it is possible to overcome the dynamic range 
limitation previously reported. The accuracy of the results given by the sMQ method for large 
distances and long times decreases with increasing absorption. For \i a = 0.001 mm -1 if we ex- 
clude the first temporal bin, the relative error does not exceed 10% across all the values of the 
reduced scattering coefficient examined. For increasing ji a values, the relative error increases 
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Fig. 6. (a) Reflectance signal obtained with EFD from the cMC simulation, (b) relative 
error of the reflectance provided by the sMC^ method, (c) relative error of the reflectance 
provided by the sMQ method, for ji a = 0.1 mm -1 and \i' s = 1.5 mm -1 . 

for large distances and late times. Larger ji' s values increase the temporal range over which 
the error remains low, but decrease the spatial range where the model is accurate. This may 
be due to the fact that a larger scattering coefficient increases the number of collisions that a 
photon undergoes before exiting the tissue thus it increases the signal at later times, while it 
concentrates the light distribution over a smaller radial range. 

As reported by Kienle and Patterson, interpolation errors are largest for short times and small 
distances and well as for ji' s values far removed from the reference ji' s value. This arises due 
to shortcomings of discretization methods that tend to be inaccurate in regions where rapid 
spatial/temporal variations of the reflectance signal exist. Figure 7 compares the the relative 
error measured for times <0.1 ns and source detector locations <15 mm for \i a = 0.1 mm -1 and 
ll' s = 0.5 mm -1 using (a) the NURBS surface built upon the reference nominal data or (b) linear 
interpolation of the reflectance binned in uniform bins with At = 5 ps and Ar = 0.2 mm. This 
result shows that the NURBS sMQ approach typically reduces the relative error by an order of 
magnitude and provides accuracy comparable to that obtained using the sMC^ approach. This 
provides a reflectance signal with excellent accuracy even when the optical properties chosen 
represent a large perturbation from the reference reflectance. 

5. Conclusion 

We have demonstrated how the scaling relations used for the single Monte Carlo (sMC) method 
can be derived from the radiative transport equation. While the scaling associated with the scat- 
tering coefficient is rigorous, the application of Beer's law to model the effect of absorption is 
only approximate when applied to binned data. As a result, the accuracy of the sMC approach 
degrades for larger ji a values and/or when the discretization intervals are large. We show that 
when sMC is applied on a photon by photon basis it yields reflectance estimates that are statis- 
tically equivalent to those obtained using conventional MC simulations. 

However a more practical implementation of the sMC requires the development of improved 
spatial and temporal binning and interpolation approaches. The challenge in such methods is 
to construct a reference reflectance in a fashion that minimizes the discretization error and 
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Fig. 7. (a) Logio of the relative error of the sMQ approach using NURBS interpolation 
on the reference values, (b) logio of the relative error of the sMQ approach based on 
linear interpolation of the uniformly binned reflectance as done by Kienle and Patterson for 
jl a =0.1 mm -1 and ji' s = 0.5 mm -1 

error propagation due to interpolation. We demonstrate the advantages of the equal frequency 
discretization approach as compared to equal width discretization, which provides a low, and 
nearly uniform, relative error over the entire measurement domain. Even better performance 
can be realized using a novel adaptive binning strategy that accommodates the large dynamic 
range of the reflectance values while also reducing interpolation errors. The regularization of the 
number of photons collected in each bin leads to an almost uniform relative standard deviation 
across all the measured reflectance values over a wide physical range. Moreover the windowing 
mechanism implemented in the adaptive strategy provides a more dense network of nominal 
values that can be robustly interpolated. 

Our implementation of NURBS interpolation provides a fast and stable algorithm that re- 
quires a small amount of information as compared to other fitting methods. The results shown 
here indicate that the sMQ method can be used to derive R(r,t) for a wide range of optical 
properties and source detector separations. Given its speed, the sMQ method can potentially 
be used as an inverse solver engine for the recovery of optical properties from experimental 
data as it provides transport-rigorous estimates, is based on very fast and stable algorithms and 
requires little data storage. The sMQ forward solver based on adaptive binning and NURBS 
interpolation (NURBS Forward Solver) is available for use in the Virtual Tissue Simulator at 
http : / /www . virtualphotonics . org/ . 
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